{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monty Hall 问题\n",
    "参赛者面前有三扇关闭的门，其中一扇门后有一辆汽车，而另外两扇门后面各藏有一只山羊。参赛者从三扇门中随机选取一扇，门后物品即奖励。当参赛者选定一扇门尚未开启时，主持人会打开一扇背后有羊的门，然后问参赛者是否要换一扇门。问：参赛者此时换门，是否会影响他赢得汽车的概率？[1]\n",
    "\n",
    "## 实验设计\n",
    "显然，不换门的胜率为 $\\frac{1}{3}$，因此我们只需模拟并统计换门的胜率，再通过比较和分析，就可以得出结论。为此，我们将首先进行 $N$ 次重复测试，然后：\n",
    "1. 统计当 $N$ 线性增长时，频率和标准差的变化；\n",
    "2. 固定 $N$，进行多组实验，统计胜率的随机分布情况，观察其是否呈正态分布。\n",
    "\n",
    "## 设计思路\n",
    "在模拟过程中会出现 $0, 1, 2$ 三个下标，已选其中一个或两个，比如一个是 `open`，代表主持人开的门，一个是 `choice`，代表参赛者选的门，如何在不读取 `open` 和 `choice`，也不遍历 $0, 1, 2$ 的情况下直接选出没有被选择的那一个数是一个能提高效率的做法。事实上，注意到不论 `open` 和 `choice` 为何值，必有\n",
    "$$\n",
    "\\mbox{open} + \\mbox{choice} + x = 0 + 1 + 2 = 3 \\Rightarrow x = 3 - \\mbox{open} - \\mbox{choice}. \n",
    "$$\n",
    "\n",
    "## 模拟结构\n",
    "**one_test()** 一次试验，返回成功与否（True/False）；\n",
    "**batch_test(N)** $N$ 次重复试验，返回成功率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def one_test():\n",
    "    c = np.random.randint(0, 3)   # 随机设置一个车的门号\n",
    "    choice = np.random.randint(0, 3)   # 参赛者选一个门号\n",
    "    if choice == c:   # 如果参赛者选的门后有车\n",
    "        open = np.random.randint(0, 2)   # 主持人有两个门可选，随机选\n",
    "        if choice <= open:   # 后移去掉参赛者已选择的门\n",
    "            open = open + 1\n",
    "    else:   # 参赛者选的是羊门，此时主持人的选择是唯一的             \n",
    "        open = 3 - choice - c\n",
    "    choice = 3 - choice - open   # 参赛者换门，非 choice 非 open\n",
    "    if (choice == c):\n",
    "        return True\n",
    "    else:\n",
    "        return False   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_test(N):\n",
    "    success = 0\n",
    "    for i in range(N):\n",
    "        if (one_test()):\n",
    "            success += 1\n",
    "    return success / N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验一：统计当 $N$ 线性增长时，频率的变化\n",
    "重复次数初始为 $100$，每次增加 $100$，一共进行 $100$ 组试验，观察频率的变化。从结果图像中可以看到，随着重复次数 $N$ 的增加，获胜的频率稳定在 $\\frac{2}{3}$ 附近，振荡幅度逐渐收窄。这一现象和我们的理论估计一致。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 100    # 总试验组数\n",
    "N0 = 100     # 试验起始规模\n",
    "cases=np.arange(N0, N0 + T * 100 + 1, 100)        # 产生全部试验规模\n",
    "sr=[batch_test(i) for i in cases]   # 逐个试验并记录"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(cases, sr, 'b.')\n",
    "plt.plot([N0, N0 + T * 100 + 1], [2/3, 2/3], 'r-')   # 参考线 2/3\n",
    "axs = plt.gca()\n",
    "plt.axis([N0, N0 + T * 100 + 1, 0.4, 0.8])\n",
    "axs.set_xlabel(\"Number of repeat\")\n",
    "axs.set_ylabel(\"Success rate\")\n",
    "plt.savefig(\"success_rate.pdf\")\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验二：固定 $N$，进行多组实验，统计胜率的随机分布情况，观察其是否呈正态分布。\n",
    "重复次数固定为 $500$，一共进行 $100000$ 组试验。绘制模拟结果的分布。从结果图像中可以看到，该分布符合期望为 $\\frac{2}{3}$ 的正态分布。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 500\n",
    "T = np.arange(100000)\n",
    "sr=[batch_test(N) for i in T]\n",
    "plt.hist(sr, bins = 16)\n",
    "axs = plt.gca()\n",
    "axs.set_xlabel(\"Success rate\")\n",
    "axs.set_ylabel(\"Number of times\")\n",
    "plt.savefig(\"hist.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 结论和进一步工作\n",
    "模拟的结果在关于重复次数增长率和大样本（每组 $500$ 次重复，共 $100000$ 组）的结果分布上，均服从期望为 $\\frac{2}{3}$ 的正态分布，从而验证了换门之后，胜率为 $\\frac{2}{3}$ 的理论结果。\n",
    "\n",
    "可以考虑的进一步工作有：\n",
    "1. 如果进步一计算标准差等，可以在二阶或更高阶矩上，验证试验结果和分布与理论结果的一致性；\n",
    "2. 可以进一步根据模拟结果，计算胜率为 $\\frac{2}{3}$ 的点估计和区间估计。\n",
    "以上问题留待后续问题解决。\n",
    "\n",
    "## 参考文献\n",
    "[1]肖柳青, 周石鹏. 随机模拟方法与应用. 北京大学出版社, 2014."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
